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Prior to recombination photons, electrons, and atomic nuclei rapidly scattered and behaved, 
almost, like a single tightly-coupled photon-baryon plasma. We investigate here the accuracy of 
the tight-coupling approximation commonly used to numerically evolve the baryon and photon 
perturbation equations at early times. By solving the exact perturbations equations with a stiff 
solver starting deep in the radiation-dominated epoch we find the level of inaccuracy introduced 
by resorting to the standard first-order tight-coupling approximation. We develop a new second- 
order approximation in the inverse Thomson opacity expansion and show that it closely tracks the 
full solution, at essentially no extra numerical cost. We find the bias on estimates of cosmological 
parameters introduced by the first-order approximation is, for most parameters, negligible. Finally, 
we show that our second-order approximation can be used to reduce the time needed to compute 
cosmic microwave background angular spectra by as much as ~ 17%. 

PACS numbers: 98.80.-k,98.80. Jk 



I. INTRODUCTION 

The cosmic microwave background (CMB) radiation 
provides us with a picture of the Universe as it looked 
when the first atoms formed, about 380 000 years after 
the big bang. At that time, photons and baryonic 
matter practically ceased interacting and the Universe 
became transparent to radiation, allowing CMB photons 
to free-stream through space. To extract accurate 
cosmological information from CMB data it is crucial 
to understand the evolution of the photon-baryon 
plasma before decoupling. This involves solving the 
Boltzmann equations for both photons and baryons 
coupled by a Thomson-scattering collision term [IHZ]- 
However, the large value of the Thomson opacity (r~^) 
before recombination renders these equations stiff, and 
hence difficult to solve numerically. This difficulty is 
usually circumvented by making use of the so-called 
"tight-coupling" approximation P^. In this scheme, 
an alternative (approximate) form of the equations is 
found and used to find the solution by systematically 
expanding the problematic terms to first order in t^. 
At late times, once the Thomson opacity drops below 
a certain threshold, one switches back to the exact 
equations to determine the final answer. 

Recently, it has been shown that uncertainties in the 
cosmological recombination process may lead to a bias in 
estimates of cosmological parameters [SHIO]. Could the 
tight-coupling approximation also result in such a bias 
and affect the final result of modern Boltzmann codes 
such as CAMB [11 or CMBFAST [12]? In this paper, we first 
investigate the accuracy of the tight-coupling approxima- 
tion by directly solving the exact set of equations at all 
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times using a stiff integration scheme. This necessitates 
calculating more accurate cosmological initial conditions 
than has been done in the past. While not efficient, solv- 
ing the exact equations allows us to determine the level of 
inaccuracy introduced by resorting to the tightly-coupled 
limit at early times. We then design a higher-order ex- 
pansion scheme and show that at second order in kTc and 
Tc, the final solution very closely tracks that obtained by 
solving the exact set of equations. We are then able to 
compute the bias on cosmological parameter estimates 
introduced by resorting to the first-order tight-coupling 
approximation and show that it is indeed small for most 
cosmological parameters. Finally, and most importantly, 
we describe how our second-order expansion can be used 
to speed up the computation of CMB power spectra with- 
out loss of overall accuracy. 



II. SOLUTION TO THE EXACT EQUATIONS 

The first step in testing the validity of the tight- 
coupling approximation is to evolve the exact set of equa- 
tions from early times. This requires the use of a dif- 
ferential equation solver able to solve stiff systems with 
adaptive step sizes. We utilize the LSODA fT^ solver 
which is based on the backward differentiation formula 
method. We find that the stiff integrator can solve the 
exact Boltzmann equations provided suitably accurate 
initial conditions are given. Indeed, the usual initial con- 
ditions for the perturbation variables used by modern 
Boltzmann codes are valid only in the limit of perfect 
coupling between photons and baryons [3 |T3] . In this 
limit, the dipole moments of the photon and baryon dis- 
tributions are exactly equal to each other and the pho- 
ton quadrupolc moment vanishes. However, in order to 
solve the exact equations at early times, one needs to 
initialize the relative dipole moment (usually called the 
slip) between the photons and baryons and the photon 
quadrupole moment to nonzero values. We describe our 
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approach to this problem in the next subsection. We 
then verify the convergence of the solution obtained with 
the stiff integrator to ensure it is stable to changes in the 
numerical tolerance and accuracy settings. 

A. Initial Conditions 

To find suitable initial conditions to the system of ex- 
act equations, we expand each perturbation variable in 
powers of kr and e = Tc/t 

A(T,e) = ^(CA)™„(fcTre" (1) 

m.n 




and substitute the result in the system of coupled dif- 
ferential equations (see Appendix |A]). Here k is the 
Fourier wave number, r is conformal time and A(r, e) 
stands for any of the following perturbation variables: 
5c,S^,0-y, F^2, Sb, Sb = Ob — 0j,5v,0u, -Fl/2, V (our notation 
closely follows that of [7]). We then match coefhcients 
of like powers of kr and e to obtain a set of linear equa- 
tions for the series coefficients (CA)mn- We then solve 
these linear equations to find a global series solution, de- 
manding that the tightly-coupled solutions (adiabatic or 
isocurvature) are retrieved in the limit e — > 0. In prin- 
ciple, one could try to solve the full recursion relation 
and obtain a closed-form expression for the (CA)mri- In 
practice however, finding the first few terms of the series 
is sufficient to set accurate initial conditions. Using this 
method, we obtain the leading-order contribution to the 
initial value of the slip between baryons and photons for 
the adiabatic mode 



Sb{T) = Obir) - e^{T) 



6(1 - 



-L0k\^€ + O{€^), (2) 



where Pi — 1 — l{l + 2)K/k'^ is a normalization constant, 

Ru = Pvl{pv + P7), ^6 = PhlPm and W = iJo^^m/V^- 
The leading-order contribution of the photon quadrupole 
moment of the adiabatic mode is 
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We list the initial conditions for all of the relevant per- 
turbation variables in Appendix [B] 



FIG. 1: Fractional change in Cf^ and Cf^ versus multipole 
moments as the relative tolerance of the stiff integrator is 
varied from lO"*^ to 10"'^. Here the GAME accuracy boost 
factors are set equal to 5. The average change is 3.5 x lO"'^ 
for Cj'^ and 3.3 x 10"'^ for Cf ^. 



hierarchies and the sampling of the final angular power 
spectrum. See Ref. [T5] for a complete list. Here, we 
increase the accuracy boost factors to verify for conver- 
gence but we also vary independently the tolerance of 
the stiff integrator to single out any error that is intro- 
duced by the solver itself. Throughout this section, we 
use as a benchmark model the WMAP seven-year cos- 
mological parameter best- values [TB]. Figure [l] shows the 
fractional change in both Cj'^ and Cf^ as a function of 
the multipole moment I as the relative tolerance of the in- 
tegrator is increased by an order of magnitude from 10~^ 
to 10"'^. The average fractional change in the angular 
power spectra is approximately 3 x 10"'' hence show- 
ing that the integration process has converged. Figure 

t shows the fractional change in both Cj^ and Cf^ as 
e three CAMB accuracy boost factors are increased from 
5 to 6. We see that the Ci computed with the stiff in- 
tegrator have an accuracy of 0.01% or better with the 
accuracy boost factors set to 5. We shall use this spec- 
trum as our benchmark for testing the accuracy of our 
second-order tight-coupling approximation scheme which 
we now present in the next section. 



B. Convergence of the Stiff Integration 

We verify the convergence of the stiff integrator by 
running several computations with increasing accuracy 
and comparing the resulting angular power spectra. In 
CAMB, the desired accuracy is usually selected by choosing 
the appropriate "accuracy boost factors" which control, 
among other things, the Fourier mode sampling of the 
CMB anisotropy sources, the time step of the integrator, 
the number of multipoles kept in photon and neutrino 



III. SECOND-ORDER SCHEME 

In the usual tight-coupling approximation, the photon 
and baryon dipole moments are obtained by solving the 
two exact equations |17] 
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FIG. 2: Fractional change in Cf^ and Cf"^ versus multipole 
moments as the three CAMB accuracy boost factors are in- 
creased from 5 to 6. The maximum change is about 1 x 10"* 
for Cf ^ and 6 x 10"^ for Cf^ . 
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where a stands for the scale factor and dots represent 
derivatives with respect to conformal time. The exact- 
ness of the solution to the above equations of motion 
depends strongly on the accuracy at which we can de- 
termine both Sb and F~^2- Current CMB Boltzmann 
codes use a first-order expansion in Tc to approximate 
the photon-baryon slip and the photon quadrupole mo- 
ment. Here, we propose a method to obtain the second- 
order corrections in to these quantities. See |18) for a 
related expansion in the context of magnetogenesis. 



A. Photon-Baryon Slip 

Our starting point is the e xact equation for the slip ob- 
tained from combining Eq. ( A8 1 and the time derivative 
of Eq. (IaTTI) 17J: 
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Usually, one sets Sb — F^2 — F^2 — and neglect the 
prefactor on the right-hand side of Eq ^ since they 
contribute terms of order and higher. However, to 
obtain an equation for the photon-baryon slip valid at 
second order in Tc, an approximation for Sb, Fj2 and 
Fj2 accurate to first order in is necessary. 



The second derivative of the photon-baryon slip is com- 
puted by taking the time derivative of the right-hand 
side of Eq. Here, we neglect terms proportional to 
(fiSb/dr^ and F^2 ■ We then use the time derivative of 
Eqs. \A1\ and (AlOl to eliminate the second derivatives 
of 5^ and Sh- h is eliminated by using the i-i component 
of the perturbed Einstein equation 



h + 2-(2fca - 677) - 2k^/3n] = -SwGa'^ ^ 3piW^6„ (7) 



- 6ri)/ 2k is the shear. Wc further eliminate 
(A8) and set F^2 — -F72 = since they 



where a = {h- 
9j using Eq. 

contribute terms of order to Sb- We finally substitute 
the time-evolution equations for the parameters R, c^, Tc 
and d/a = H: 



R= -n{l-icl)R, 



= -net 



n = -mn - w 



Now armed with an expression for Sb, we substitute it 
into Eq. ^ and solve algebraically for Sb- The result is 
given in Appendix [C} 



B. Photon Quadrupole Moment 

To obtain an expression for F^2 and F^2 accurate to 
second order in t^, we use the recursion relation between 
higher photon multipole moments [J 



F. 



[lF,(i-i) -{1 + mF^n+i)] - ^F,i, (8) 



which is valid for I > 3. We begin by setting F^^ — 
and solve Eq. (|8| with I — 4 lor F^^. We then take 
the derivative with respect to proper time, setting F^^ — 
0. We finally solve the resulting equation for F^4 and 
substitute back the result in Eq. ([s]). This last equation 
leads to an expression for Fj4 valid to fourth order in Tc 
(remembering that _Fly3 cx fc^r^ and that Tc oc Tc/t): 



F,4 ^ \kTcF,:,{l - Tc) - ^kT^F^, + 0{tI)- (9) 



Substituting the above in Eq. ([8| evaluated at Z = 3 and 
using a similar procedure, we obtain an expression for 
F^j valid to fourth order in 

i^73 - ^kTcF^2{l-Tc + fl)-hTlF^2{l-Tc) 
1 fi 

- — k^TlF,2+0{Tl). (10) 

The last step in deriving expansions for the quadrupole 
moment and its derivative is to express the polarization 
multipole G^2 in terms of F~f2 and F^2- Similar to the 
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above calculation, this is accomplished by using the re- 
cursion relation for the polarization multipole moments 
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where Sij is the Kronecker delta. Again, we set G^5 = 
and follow the method outlined above to obtain 
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-^fcV2F^2(l - 6f,) + i^fcV^F^s, (12) 

which is accurate to fourth order in Tc- We now have all 
the necessary tools to derive approxi ma te exp ress ions for 
Fy2 and Fj2- We substitute Eqs. (12) and (10 1 in Eq. 
(A9) and solve for F^2- We then takethe derivative with 



respect to conformal time and set Fj2 — 0. We finally 



solve for F^2 and obtain 
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Substituting the above back in Eq. (A9) we ulti- 



mately arrive at the desired expression for the photon 
quadrupole moment 
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C. Computational Procedure 



As we can see from Eq. ( 14 ), our second-order expres- 



sion for the photon quadrupole moment depends on 9^. 
From a practical perspective, this is problematic since it 
is the quantity that we are trying to determine in the 
first place. We overcome this difficulty by computing 
each quantity order by order in until the desired level 
of accuracy is reached. 

The first step is to obtain an approximation to F^2 
using Eq. (14 1 but keeping only the terms linear in r^, 
which are independent of 9^ . We then use this expression 
to calculate & to first order in using the traceless space- 
space part of the perturbed Einstein equation 

k& + 2°-ka - k^T] = -SttGo^ {pi,F^2 + P',F^2) ■ (15) 

Next, we calculate a zeroth-order expression for 9j using 
Eq. (U| with Sh and Fj2 set to zero. We then use our 
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order tight-coupling approximation and the benchmark in- 
tegration (full line), and between the second-order approxi- 
mation and the benchmark integration (dashed line). Here, 
the three sets of C;s have been computed with default accu- 
racy. The average fractional difference is 6.6 x 10-* for the 
first-order approximation and 5 x 10~^ for the second-order 
approximation. 



two formulas for a and 9-y to compute F^2 to first order 



in Tc using Eq. (13). 



We now have all the necessary tools to calculat e th e 
photon-baryon slip to second order in Tc using Eq. (CI ). 



We finally use this last express ion t o obtain a first order 
approximation to 9^ using Eq. ( A8 ) which in turn is used 



to obtain Fj2 and Fj2 accurate to second order in Tc 



D. Accuracy of the Second-Order Scheme 

We test the accuracy of the second-order scheme by 
comparing the final angular power spectrum with both 
the stiff integrator benchmark and the usual first-order 
tight-coupling approximation. To isolate the effect of 
the second-order terms in the equations of motion, we 
leave untouched the algorithm that switches from the 
tightly-coupled equations to the exact system of equa- 
tions. Improvements to the switching criteria will be dis- 
cussed in Sec. |V] As mentioned above, all the results 
presented in this section are valid for the WMAP seven- 
year best-fit values for cosmological parameters. We find 
that at default accuracy setting ("accuracy boost" — 1) 
for all three computations, the fractional difference be- 
tween the second-order scheme and the benchmark in- 
tegration averaged over multipoles is about an order of 
magnitude smaller than the average fractional difference 
between the usual first-order tight-coupling approxima- 
tion and the benchmark integration (see Fig. [3|. Hence, 
the second-order scheme reproduces more accurately the 
solution to the exact equations. 

As the accuracy boost factors are increased, the 
second-order scheme keeps providing, on average, a more 
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FIG. 4: Fractional difference of Cf^ between the usual first- 
order tight-coupling approximation and the benchmark inte- 
gration (full line) and the second-order approximation and 
the benchmark integration (dashed line). Here, all the Cis 
have been computed with the three accuracy boost factors 
set to 5. The average fractional difference is 2.4 x 10~^ for 
the first-order approximation and 3.5 x 10~^ for the second- 
order approximation. 



accurate answer than the first-order tight-coupling ap- 
proximation. Figure |4] compares the angular power spec- 
tra from the two schemes with those found by integrating 
the exact system of equations. Although, the difference 
between the two codes might be insignificant for current 
CMB experiments, it illustrates that the next-to-leading- 
order code is better capturing what is happening during 
the tightly-coupled epoch, especially for the low multi- 
poles. The key point however is that this better accu- 
racy comes at almost no additional computational cost, 
a point that we shall emphasize in Sec. [V[ 

In summary, we have shown that the second-order 
tight-coupling approximation reproduces more closely 
the result found by solving the exact equations, hence 
showing that the tight-coupling expansion is converging 
toward the exact solution. For practical applications 
however, the percentage change in the angular power 
spectrum between the usual first-order approximation 
and the exact solution is small and well within the quoted 
precision from CAMB (0.3% at default accuracy). This 
implies that the first-order tight-coupling approximation 
should be sufficient for most practical purposes. Never- 
theless, as we will describe in the next few sections, it is 
possible to use our second-order expansion to reduce the 
potential bias on cosmological parameter estimates and 
to speed up the code. 

We mention in passing that the precision (i.e., the size 
of the numerical noise) of individual Ci is almost not af- 
fected by the introduction of the second-order terms in 
the equations of motion. Indeed, the precision of the fi- 
nal angular power spectrum is mainly determined by the 
number of k-modes evolved by the code, the number of 
photon multipoles that are solved for, as well as various 
interpolation errors. Since our new tight-coupling ap- 



proximation does not modify any of the above, it is there- 
fore natural to expect that the precision of the second- 
order power spectrum to remain unchanged. 



IV. BIAS ON COSMOLOGICAL PARAMETERS 

In today's era of precision cosmology, the main purpose 
of CMB codes is to generate theoretical spectra that are 
then compared with data for cosmological parameter es- 
timation purpose. However, numerical errors in the the- 
oretical spectra could lead to a slight bias on estimates of 
cosmological parameters jl5] . Since our improved tight- 
coupling approximation scheme leads to more accurate 
values of the power spectra, it is interesting to see how 
the bias is affected. To answer this question, we need to 
compute the effective between our theoretical spectra 
and a fiducial data set which we take to have Planck-level 
noise. The effective is defined by 
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where /sky is the observed sky fraction and Ci = {C( 
M{^'^ } is the theoretical covariance matrix. Here, X 
runs over temperature (T) and polarization (E). Ci is 
the data covariance matrix. If we assume that the like- 
lihood C = exp {—x^/2) is a multivariate Gaussian and 
that the prior probability densities are flat, then the bias 
on any cosmological parameter cannot exceed i/x^ stan- 
dard deviations. In practice however, this bound is rarely 
saturated. Nonetheless, a small x^ between the data and 
the theory is still necessary to ensure a minimal bias. 

We generate a fiducial data set up to / = 2500 using the 
method outlined in p[5j but with the C; obtained from 
the stiff solver. Again, we use the WMAP 7-year best-fit 
values for cosmological parameters. We take the noise to 
be Gaussian and isotropic with power spectrum given by 
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where ^beam is the beam width and Ax is the sensitivity 
per pixel. As an example, we consider the 143 GHz chan- 
nel of the HFI instrument aboard Planck |^ which has 
dhcam = 7.1', At — G-OfiK and A^ = llAfiK, assuming 
14 months of integrated observation. We assume a sky 
coverage of /sky = 0.65. 

We list in Table |l] the values of x^ computed between 
our fiducial Planck data and our improved second-order 
code. For comparison, we also give the x^ values for 
unmodified CAMB at similar accuracy boost. We see that 
the higher-order tight- coupling approximation leads to a 
better fit to the fiducial data and therefore to a smaller 
theoretical maximal bias on cosmological parameters at 
no extra numerical cost. To estimate the real biases on 
cosmological parameters, we run several Markov chains 
using both the first- and second-order tight-coupling code 
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together with the pubHcly available code CosmoMC TD^. 
We restrict ourselves to the "vanilla" G-parameter ACDM 
model and made sure that the Gelman-Rubin "R — 1" 
convergence criteria [5^ is smaller than 0.005 for all the 
parameters under consideration. 

We list in Table |IT] the biases between the results from 
our second-order CMB code and the results from a code 
that used the same accuracy setting as the fiducial spec- 
tra (mimicking an error-free analysis). For comparison, 
we also give the biases for the usual first-order code. At 
default accuracy, we see that the difference between the 
two codes in terms of the bias is rather slim, with 9 being 
the most dramatically affected by the second-order code. 
This stems from the fact that our second-order code bet- 
ter captures the position of the first peak as can be seen 
from Figs. [3]and|4] We conclude that the numerical er- 
rors due to the first-order tight-coupling approximation 
does introduce a small bias to the estimate of 6 at default 
accuracy, although it is clear that other numerical errors 
(k-sampling, interpolation, etc) contribute the most sig- 
nificant part to the biases of cosmological parameters for 
both codes. As the accuracy of the codes is increased, 
the difference in bias between the two codes becomes in- 
significant for parameter estimation purposes. Therefore, 
if one sets the accuracy of the theoretical spectra to be 
large enough such that the bias from numerical errors 
others then the first-order tight-coupling is small, then 
the usual tightly-coupled equations are appropriate for 
cosmological parameter estimation. 



V. REDUCING THE COMPUTATIONAL 
RUNTIME 

Up to this point, we have used the second-order ex- 
pansion in Tc to improve the accuracy of CMB Boltz- 
mann codes. In this section, we adopt a different point of 
view and take advantage of our improved tight-coupling 
scheme to reduce the computational time needed to 
evolve the perturbation equations. Indeed, the new 



Code 


x' 


Time (s) 


CAMB accuracy — 1 


2.3 


4.8 


2nd Order acc. = 1 


1.3 


4.8 


CAMB accuracy — 2 


0.17 


18.2 


2nd Order acc. = 2 


0.091 


18.2 


Opt. CAMB acc. =2 


1.1 


15.1 


Opt. 2nd Order acc. = 2 


0.10 


15.1 



TABLE I: values between fiducial Planck data and theo- 
retical spectra gotten with the first- and second-order codes 
for different accuracy boost. We also give the computational 
time needed to generate the theoretical spectra in order to 
show that the greater accuracy comes at no extra numeri- 
cal cost. The computational times displayed here are for a 
single-processor machine. 



Code 






e 


r 




ln(10"M,) 


CAMB accuracy — 1 


0.24 


0.15 


0.31 


0.11 


0.40 


0.19 


2nd Order acc. = 1 


0.25 


0.16 


0.15 


0.12 


0.33 


0.15 


CAMB accuracy — 2 


0.02 


0.006 


0.03 


0.013 


0.03 


0.03 


2nd Order acc. = 2 


0.03 


0.003 


0.01 


0.016 


0.017 


0.015 



TABLE II: Biases of the 6-parameter ACDM model in unit 
of the standard deviation. We contrast the first- and second- 
order tight-coupling approximation and give the value of the 
accuracy boost factors used for each computation. 



0(t^) terms in the tightly-coupled perturbation equa- 
tions allow one to switch to the exact equations at a later 
time while keeping the same accuracy as the usual first- 
order expansion. Since the approximate tightly-coupled 
equations are easier to solve than their exact counter- 
parts, precious computational time can be saved. More- 
over, the higher accuracy of the second-order equations 
lets us use a larger minimal time step for the numerical in- 
tegrator, hence reducing the total number of steps taken 
by the integrator and further cutting down the compu- 
tational time. 



Our approach here is to degrade the accuracy of 
the second-order code by modifying the tight-coupling 
switching criteria, using larger time steps and cutting 
down the photon hierarchy until the output from this 
"optimized" code somewhat matches that of the unmod- 
ified first-order code. We then compute the value be- 
tween our fiducial Planck data and the output from this 
optimized code and compare it to a similar calculation 
done with regular CAMB. The results are shown on the 
two last rows of Table |T] where we see that we achieve a 
~ 17% computational time reduction while still retaining 
the accuracy of the first-order code at accuracy boost 2. 

Although this gain in computational efficiency is mod- 
est, it can significantly reduce the amount of time neces- 
sary to run Markov chains for cosmological parameter es- 
timation. To demonstrate this, we run 8 chains with our 
optimized second-order code at accuracy boost 2, gener- 
ating 20000 samples per chains. We also run the similar 
chains with regular CAMB at accuracy boost 2 and make 
sure to have -R — 1 < 0.005. Figure [5] shows that the re- 
sults for the marginalized posterior distribution are very 
similar between the two codes, with the distribution for 9 
being the most affected, although very mildly (0.09 stan- 
dard deviation). However, the most important difference 
between the two results is that our optimized second- 
order code took on average ^ 16% less time to complete. 
Hence, our second-order tight-coupling code, in addition 
to leading to more accurate angular CMB spectra, can 
instead be used to speed up the computational time and 
make more efficient use of computing resources. 
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FIG. 5: Marginalized posterior probability distribution for 
the vanilla ACDM model. The full black line represents the 
result gotten using the first-order code at accuracy boost 2 
while the red dotted line was obtained using our optimized 
second-order code. 



In conclusion, we have investigated the accuracy of the 
tight-coupling approximation by solving the exact equa- 
tions at all times using a stiff numerical solver. We have 
shown that the first-order tight-coupling approximation 
leads to a small accuracy lost compared to the exact so- 
lution and that this change is well within the quoted pre- 
cision of the angular power spectra. We have discussed 
how our second-order code has a smaller maximal possi- 
ble bias on cosmological parameters than its first-order 
counterpart. We have shown that the bias introduced by 
the first-order tight-coupling is insignificant for today's 
cosmological experiments, unless CAMB's default accuracy 
is used. Finally, we have shown that the improved ac- 
curacy of our second-order approximation allows one to 
optimize the tight-coupling switching criteria and inte- 
gration parameters in order to reduce the computational 
time of the code. (After this project was completed, a 
related work appeared [23].) 
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FIG. 6: Residuals (A(5^ = S^""""' - between the pho- 

ton perturbation S-f computed using the exact equations and 
the solutions obtained with the first- and second-order tight- 
coupling approximation. 



VI. DISCUSSION AND CONCLUSION 

We have developed a second-order tight-coupling ap- 
proximation to the photon-baryon perturbation equa- 
tions and shown that it closely reproduces the solu- 
tion to the nonapproximated equations. In practice, the 
main reason why our second-order tight-coupling code 
produces more accurate power spectra is that it better 
tracts the evolution of the photon perturbations. Figure 
[6] shows the residuals between the photon perturbation 
6j computed using the exact equations and the solutions 
obtained with the first- and second-order tight-coupling 
approximation. We see that the second-order scheme de- 
viates much less from the exact solution then the usual 
first-order scheme, leading to a more accurate value of 
the source term needed for the line-of-sight integration 
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Appendix A: Perturbation Equations 

In this Appendix, we list the perturbation equations 
used to solve for the initial conditions found in Appendix 
[b| We closely follow the notation of [7]. Here, r] and h 
stand for the synchronous gauge curvature perturbation 
variables, a is the scale factor, K is the inverse of the 
squared curvature radius, pi is the energy density of the 
j"* specie, Wi = Pi/ pi, where pi is the pressure, and a dot 
denotes differentiation with respect to conformal time. 

k'^f3iV~l-h + 4TTGa^y" pA = (Al) 

i 

Kh 

k'Pif]- — -^^Ga^Y.^l + w,)pA = (A2) 

i 

5c + \h = Q (A3) 
5u + \o, + \h = (A4) 

- ^ {5u - 2/3iF,2) = (A5) 
- + ^/32fcF,3 - ^ ( J + J =0 (A6) 
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h 
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4 Tf, 



(A7) 
(A8) 
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15 V 4 
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72 



G 



72 



Sb — 



\ + R 



(Sb + O-y) — 5*6 



a 
1 



■A 



72 



(A9) 
O(AIO) 

O(All) 



Here, R = (4/3)p^/ pb- Our approach to solve these equa- 
tions follow closely that of [24^. We first use Eq. (Al) 
to eliminate h in favor of the curvature perturbation rj. 
For simplicity, we set = 0. We then approximate the 
octupole moment of the neutrinos and photons as: 



F,, 



kr 



315 



(A12) 



1. Photons 



F. 



73 



-kTcF^2- 



(A13) 



Finally, we eliminate the photon polarization moments 



from (A9) using 
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F. 



5tc 
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F. 



72 
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Appendix B: Initial Conditions 



In this Appendix, we list the initial conditions ob- 
tained by the method outlined in Sec. |II A[ Here, 

Ru = Pu/{Pu + P-t), Rb = Pb/Pm, W = Hofl,n/V^, 

e — Tc/t, Sbir) = Ob{T) — 9j{t) is the velocity difference 
between baryons and photons. Note that our convention 
for the normalization of perturbations differs from |24| 
by /3i — —(3i/2. Note also that what we label Pi here is 
denoted by (32 in [24) . 



S,ir) = - — k r + - 



3 + A(4Ai^. + m-5) 
27 (4i?, + 15) 



/^l 2 I 2 4 

— w k T 
24 



(Bl) 



18 36i?„ + 135 



k\h 



\ (1 + bRb - Ry) 



120 (1 - R^ 
2(3i (2 {5Rb - 9) i?^ + 75i?h + 8Rl + 10 



15 (i?^ - 1) {2R^ + 15) (4i?^ + 15) 



I 4 4 



16/3i (6i?, + 181) 4 3 2 



45 (2i?^ + 15) (4i?„ + 15) 



(B2) 



J^72(r) 



64 



9(4i?^ + 15) 



fcVe 



4 {8R, - 5) 



3 (2i?„ + 15) (4i?^ + 15) 



32 (6i?, 181) ^2^2 ^2 



16 {2R^ {12R^ + 767) - 1855) 2 3 2 
' 9 (2i?^ + 15) (2i?^ + 25) (4i?^ + 15) ^ 



9 (2i?^ + 15) (4i?^ + 15) 



(B3) 



2. Baryons 



^bir) = -^7(t) 



(B4) 



A^,„7.4„4 
(1-i?.) 



10f3iRb 
3(l-i?,)(4i?, + 15)' 
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3. Cold Dark Matter 
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4. Neutrinos 
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5. Curvature (synchronous gauge) 



r?(r) = 2+ 



12i?j. +45 6 
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9 (2i?, + 15) (47?, + 15) 
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Appendix C: Tight-coupling Approximation to Second Order in Tc 

In this Appendix, we give the key result of our improved tight-couphng approximation scheme: the photon-baryon 
slip to second order in t^- 



re .,2 

^ nj 

Tc 1 + R 



a 

-t 
a 



1 + R 

2R (s-H^c^ + {R+ 1)% - 3-^2 j 
{R + lf 



SbTc 



{l + RY 



an{{2-3cl) R~2)eb Hk^l - 3c^)9b 



{R+1) 



+ 



3(i?+l) 



a Pc^Jb fc4(3c2 - l)c2,55 k^R{3c'i^l)S., aP{2 + ?,R)5^ U^W {{2 ^ ic^^) R - l) 5. 



a(i?+l) 3(i?+l) 



12(i?+l) a 4(i?+l) 



+ 



2(i?+l) 



^ Hfc2c^(l + (3c^-2)i?)4 ^ m'' (2 + (5 - 3c^) j?) 5^ ^ 2H(1 - 3c^)fc3a ^ fc4(3^2 _ 



+2^2(30^-1)77 + 



i? + l 4(i?+l) 
e{l - 3c2)A 



+ 



416*6 - Ak'^c% + 2Hk'^5^ + fc^^^ 
2(i?+l)2 
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Here, A = SirGa^ {p^5j + + Sc^pb^t) 
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